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ABSTRACT 



We have performed a three-dimensional magnetohydrodynamic simulation to study the emer- 
gence of a twisted magnetic flux tube from -20, 000 km of the solar convection zone to the 
corona through the photosphere and the chromosphere. The middle part of the initial tube is en- 
dowed with a density deficit to instigate a buoyant emergence. As the tube approaches the surface, 
it extends horizontally and makes a flat magnetic structure due to the photosphere ahead of the 
tube. Further emergence to the corona breaks out via the interchange-mode instability of the pho- 
tospheric fields, and eventually several magnetic domes build up above the surface. What is new 
in this three-dimensional experiment is, multiple separation events of the vertical magnetic ele- 
ments are observed in the photospheric magnetogram, and they reflect the interchange instability. 
Separated elements are found to gather at the edges of the active region. These gathered elements 
then show shearing motions. These characteristics are highly reminiscent of active region obser- 
vations. On the basis of the simulation results above, we propose a theoretical picture of the flux 
emergence and the formation of an active region that explains the observational features, such as 
multiple separations of faculae and the shearing motion. 

Key words. Magnetohydrodynamics (MHD) - Sun: chromosphere - Sun: corona - Sun: interior 
- Sun: photosphere - Sun: surface magnetism 



1. Introduction 

Solar active regi ons incl uding sunspots are generally thought to be the consequence of the magnetic 
flux emergence (Parked 1 1 95 5l) . Observationally, they appe ar as a developing bipolar pair in the 



photospheric magnetog ram and arch-filament system in Ha (lZwaanfl985l) . 



Strous et al 



(1996) and 



Strous & Zwaanl (Q999) find the hierarchy of motions of magnetic elements in the active region: 



faculae of positive and negative polarities separate from each other toward the edges of the region 
while pores of each polarity move along the edges toward the major sunspots. 

Dynamics of the flux emergence have been stu died widely through extensive numerical Simula 
tions, both in two and three dimensions (e.g. 



ShibataetalJ 1989; 



an 



Fan 2001 



Archontis et al. 2004). 
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It has been shown that the flux tube rises due to the Parker instability (IParkerl 1966) and expands 
into the corona in a self-similar way, which explains the characteristics of the flux emergence such 
as £2-shaped loops and the downflows at the footpoints of the loops. 

In this paper, we report first results of the large-scale three-dimensional magnetohydrody- 
namic (3D MHD) simulation on the emergence of a twisted flux tube from a deeper convection 
zone (-20, 000 km) into the cor ona through the photo s phere . In our previous 2D experiments, 
Toriumi & Yokoyamal {2010) and lToriumi & Yokoyamal (1201 II) (hereafter papers I and II, respec- 
tively), we found the conditions of the magnetic flux for reasonable emergence such as the field 
strength, the total flux, and the twist of the tube at a depth of -20, 000 km. We apply these val- 
ues as the initial conditions in the 3D calculation. As a result, several characteristics are found to 
be consistent with the previous 2D experiments and observations. On the basis of this simulation, 
we also suggest a new theoretical picture of the flux emergence and the birth of an active region 
through the surface. 



2. Numerical Setup 



The basic MHD equations solved in this simulation in vector form are as follows: 
| + V-(pV) = 0, 



(i) 

(2) 



dB 

dt 



V x (V x B), 



(3) 



+V- 

and 
U = 



Lu + p+±pV 2 \v+-^ExB 



1 p 



-pg-V = 0, 



y-lp 



(4) 



(5) 



1 



VxB, 



(6) 



P = —pT, 

m 



(7) 



where U is the internal energy per unit mass, / the unit tensor, &b the Boltzmann constant, m 
the mean molecular mass, and g the uniform gravitational acceleration. Other symbols have their 
usual meanings: p is for density, V velocity vector, p pressure, B magnetic field, c speed of light, 
E electric field, and T temperature. Here, the medium is assumed to be an inviscid perfect gas 
with a specific heat ratio y = 5/3. The equations are normalized by the pressure scale height 
Hq = 200 km for length, sound speed C s o = 8 km s" 1 for velocity, To = Hq/C s q = 25 s for time, 
po = 1.4 x 10" 7 g cm -3 for density, all of which are typical values for the photosphere. The units 
for pressure, temperature, and magnetic field strength are po = 9.0 x 10 4 dyn cm -2 , To = 4000 K, 
and #o = 300 G, respectively. 
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In this experiment, we set 3D Cartesian coordinates (x,y,z) 9 where z is parallel to the grav- 
itational acceleration, i.e., g = (0,0, -go)- The numerical domain is (-400,-200,-200) < 
(x/H ,y/H ,z/H ) < (400,200,250) and the total grid number is 1602 x 256 x 1024. The grid 
spacing for x-direction is Ax /Ho = 0.5 (uniform); for the other directions, Ay/Ho = 0.5 and 
Az/Ho = 0.2 at the central area of the domain, which gradually increase for each direction. We as- 
sume periodic boundaries for both horizontal directions and symmetric boundaries for the vertical 
direction. 

The background stratification consists of the adiabatically stratified convective layer (z/Hq < 
0), the cool isothermal (T(z)/Tq = 1) photo sphere/chromo sphere (0 < z/Ho < 10), and the hot 
isothermal (T(z)/Tq = 100) corona (z/Hq > 10). Here, the photosphere is isothermally stratified, 
i.e., convectively stable against the rising motion of the plasma. A transition region connects the 
low- and high-temperature atmospheres with a steep temperature gradient at around z/Hq = 10. 
The initial flux tube embedded at z = -100/Zo = -20, 000 km is uniformly twisted, i.e., B^(r) = 
qrB x (r), where B x (r) = # tu b e exp (~^ 2 /^ ube ) is the axial field, r the radial distance from the axis, 
#tube = 67£ = 2.0 x 10 4 G the field strength at the axis, R tuhQ = 5H = 1000 km the typical 
tube radius, B<p(r) the azimuthal component, and q the twist parameter which is set to be 0.1 /Ho = 
5.0 x 10" 4 km -1 (stable to the kink instability). The total flux of the axial component is <D = 6.3 x 
10 20 Mx. The middle of the tube around x/Ho = is made buoyant through thermal equilibrium 
between the inside and the outside of the tube, and the buoyancy decreases as a function of the 
form exp(-x 2 //l 2 ), where A/ Ho = 400. The tube strength, the total flux, and the twist are chosen to 
satisfy the conditions found in papers I and II. The initial background distribution of gas pressure, 
density, and temperature, and the magnetic pressure are indicated in FigureQIa). The density inside 
and outside the tube along x-axis are shown in Figure QIb). One can see that, although the density 
is reduced throughout the whole tube, the middle part is most buoyant. 

3. Results 

3.1. General Evolution 

Figure [2] shows the time evolution of the emerging twisted flux tube. Here, the logarithmic field 
strength log 10 (\B\/Bo) only in the region x/Ho < and y/Ho > is plotted. At t/ro = (Fig. 
[2ja)), the initial tube is placed at z/Ho = -100. Due to the buoyancy of the tube itself, it rises 
through the convection zone (t/ro = 400, see Fig.[2tb)). The tube expands as the external density 
decreases with height. While the plasma draining from the apex to both feet accelerates the rising 
tube, the aerodynamic drag decelerates the tube, since the external flow around the cross-section 
forms a wake behind the main tube. One can see a vortex roll at the flank of the tube and an elon- 
gated tail below. However, the azimuthal component of the flux tube yields the inward curvature 
force to maintain the tube coherent. The upper surface of the rising tube becomes fluted due to 
the interchange-mode of the magnetic buoyancy instability (magnetic Rayleigh-Taylor instability). 
Figure [3] is a two-dimensional (x-z) slice of the rising flux tube, which shows the fluting at the top 
of the tube. The rise velocity levels off while in the middle of the convection zone. 

Approaching the surface at t/ro = 600 in Fig. [3c), the tube decelerates and expands laterally 
in the v-direction to make a flat structure just beneath the photosphere (60 Mm x 20 Mm). The de- 
celeration and the flattening occur because the plasma between the rising tube and the convectively 
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Fig. 1. (a) One-dimensional (z-)distributions of the initial background gas pressure (thick solid 
line), density (dashed line), and temperature (dash-dotted line). The magnetic pressure p m = 
B 2 /(&7i) along the vertical axis x/Hq = y/Ho = is overplotted with a thin solid line, (b) 
Background density profile at z = -IOOHq (solid line). Density inside the tube along (y,z) = 
(0, -100//o) is overplotted with dashed line. 



stable photosphere is compressed, which in turn suppresses the rising tube from below. It should be 
noted that the deceleration and the flat magnetic structure are consistent with the previous 2D ex- 
periments (papers I and II). Here, the outermost field lines of the flat tube, i.e., fields at the surface 
are mostly in the negative y-direction. 

As the surficial field strength increases to s atisfy the criterion for the interchange-mode of mag- 
netic buoyancy instability (e.g. Newcomblll961 ). the secondary evolution to the upper atmosphere 
takes place at around t/ro = 660 (Fig. [3d)). We find that several magnetic domes have been built 
in the central area around x/Hq = y/Ho - at this stage, aligned in the x-direction, each be- 
ing directed in the ^-direction. As the central domes develop, another several domes are newly 
created beside the central ones (t/ro = 750; Fig. [2e)). The domes continue growing and merge 
with each other in the corona. After t/ro = 700, the rising velocity declines again, and eventually 
the emerging flux reaches z/Hq ~ 35 at t/ro - 800 (Fig.[2f)). The subsurface structure extends 
-200 < x/H < 200 and -40 < y/H < 40. 

The whole evolution process is the same as the "two-step emergence" model observed in the 2D 
experiments (papers I and II). However, the final height of the coronal structure is different from the 
2D cases: the height in the present calculation is only z/Hq = 35, while the fluxes in 2D were more 
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Fig. 2. Time-evolution of the twisted emerging flux tube. Logarithmic field strength log l0 (\B\/ Bq) 
for six different times in a limited region (x/Hq < 0,y/Ho > 0) is shown. The level of z/Hq = 
corresponds to the boundary between the convection zone and the photosphere. 




-50 
x/H 



Fig. 3. Two-dimensional (x-z) closeup of the rising flux tube at t/ro = 550. Logarithmic field 
strength log 10 (|5|/Z?o) is plotted, while the photospheric height z/Hq = is indicated with a solid 
line. 



than 200. The difference occurs due to the three-dimensionality. That is, in the 2D simulations, 
the magnetic pressure expands the magnetic flux only in the (x, z) or (y, z) plane. In the 3D case, 
however, the magnetic pressure inflates the flux in all (x, y, z) directions, whic h results in lower 
magnetic domes. This 3D situation is the same as the free expansion regime by 
Jl993L 



Matsumoto et al 
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Fig. 4. (Top) Time-evolution of vertical magnetic field strength B z /Bo at the surface z/Hq = 
(magnetogram). White (black) indicates the positive (negative) polarity. (Bottom) Corresponding 
velocity fields (Dopplergram). The vertical velocity V z /C s o (color; red is downward) and the hori- 
zontal velocity Vh/C s o (arrows) are shown. Corresponding vertical field B z /Bq at the photosphere 
is over-plotted with contours. 



3.2. Magnetic Structures in the Photosphere 

Figure |4] shows the vertical field strength (B z /Bq; magnetogram) and the corresponding vertical 
velocity (V z /C s o; Dopplergram) with the horizontal velocity field (Vh/C s o) at the surface z/Hq = 
at t/ro - 660 and 800. In the Dopplergram, red indicates a downward motion (V z /C s o < 0). At 
t/To ~ 610, we observe that some blueshifts and divergent flows appear before the flux emergence, 
which indicates that the compressed plasma between the rising tube and the photosphere escapes 
laterally at the surface. 

At t/ro = 660, in Fig.^a), magnetic elements of positive and negative polarities emerge onto 
the surface. The absolute field strength of each polarity is more than a hundred Gauss. In Fig. 
|4jc), one can find the blueshifts of the order of a few km s" 1 between each pair and the redshifts 
up to 0.12C s o = 1 kms -1 in the core of each patch. At this time, the horizontal speed in the 
positive and negative polarities is at its peak (4-8 kms -1 ), showing separative motions. These 
features indicate that the magnetic flux emerges upward, while the plasma drains downward to 
both footpoints along the field lines. Here, the surface field is mostly directed in the negative y- 
direction, and the wavelengths are A\\ ~ 20Ho and A± ~ 15 Ho respectively, where A\\ and A± are 
the wavelengths parallel and perpendicular to the surface horizontal field. The parallel wavelength 
A\\ ~ 20Hq is the most unstable wavelength of the linear Parker instability at the photosphere. 
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We speculate that the wavelength perpendicular to the field A ± is determined by the wavelength 
of the interchange-mode instability of the flux tube before it reaches the surface. During its ascent 
within the convection zone, the surface of the flux tube is found to be fluted due to the interchange 
instability (see Sect. 13.1b . When the density smoothly increases upward between the magnetized 
and the unmagnetized atmosphere, the growth rate of the interchange instability levels off as the 
wavenumber increases, and the typical wavenumber of this saturation is approximately an inverse 
of the density transition scale ( ChandrasekhaJ[l961 ) . Considering the original tube is assumed to 
have a Gaussian profile, i.e., B x (r) = # tu b e exp (-r 2 /R 2 uhe ), and thus the density transition is also a 
function of exp (-r 2 /R 2 uhQ ), the scale of this transition layer is ~ ^ tu b e . Therefore, the wavelength of 
this instability within the convection zone is approximately several times the original tube's radius 
CKtube = 5//o), which results in the wavelength of the secondary emergence A ± ~ 15 Ho at the 
photosphere. 

At t/ro = 800, the magnetic pairs develop and the region extends -120 < x/Ho < 120 and 
-30 <y/H < 30,i.e., 48 Mmxl2Mm(Fig.|4l;b)). The direction of separations is found to be tilted. 
Here, each separated patch has formed a tadpole-like shape; \B Z \ up to 350 G. This configuration 
is formed because the newly emerged elements separate outward and catch up with the elements 
that emerged earlier, and stop at the edge of the region. Also, the heads of these tadpoles make two 
alignments at the edges, and they show shearing motion. The shearing is leftward where y/Ho > 0, 
and rightward where y/Ho < 0. In this phase, the total unsigned flux JT =Q \B Z \ dxdy reaches up to 
~ 3.3 x 10 20 Mx. In Fig.|4jd) the redshifts are no longer seen, which indicates that the emergence 
has stopped. To clarify the difference between the aligned large elements (tadpole-heads) and the 
tilted elements (tadpole-tails), we use the term "pores" for the heads of the tadpole-like features. 
Moreover, this term is used in the sense of the accumulated fields that do not reach the size of a 
sunspot. 

Figure[5] shows the selected field lines plotted on the surficial magnetogram at t/ro = 800. Here, 
emerged field lines in the corona connect positive and negative polarities in the photosphere. The 
footpoints of coronal fields stop at the edge of the region, which is determined by the extension 
of the flat magnetic structure beneath the surface. That is, the size of the active region (48 Mm x 
12 Mm in this case) depends on the subpho to spheric flux tube (60 Mm x 20 Mm). The wavelength 
of the initial density deficit, A, would be one of the parameters that determine the size of the 
subphotospheric structure. We also found that some field lines connect different patches deep under 
the surface (see Figure [6]). Although such field lines are not undulating at around the surface, they 
are still reminiscent of a sea-serpent configuration and a corresponding resistive emergence process. 

The tilt of magnetic elements in the central area is caused by the emergence of the inner field 
lines (kinematic effect). Here, the initial flux tube is uniformly twisted, and thus the pitch angle 
of inner fields are smaller. Therefore, the footpoints shift as the inner field lines rise. Also, the 
shearing of two aligned "por e s" is due to the Lorentz force acting on the surface field (dynamic 



effect; see 



Manchester (2001); 



Fan 



(2001)). 



4. Summary and Discussion 

In this paper, we performed a 3D MHD simulation of the emergence of a twisted flux tube. The 
initial tube at -20, 000 km in the convection zone has a field strength of 2.0 x 10 4 G, a total flux 
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Fig. 5. Perspective view of selected field lines that pass through the vertical axis (x/Ho,y/Ho) = 
(26, 0) at t/To = 800, plotted on the photospheric vertical field B z /Bq. 




Fig. 6. A single field line that connects different photospheric magnetic elements. Plotted is the 
photospheric vertical field B z /Bq at t/ro = 800, which is transparentized. 

of 6.3 x 10 20 Mx, and a twist of 5.0 x 10" 4 km -1 , which starts rising due to its own magnetic 
buoyancy. On reaching the surface after t ~ 2.8 hr, the tube decelerates and extends horizontally 
(60 Mm x 20 Mm) owing to the convectively stable photosphere ahead of the tube. As the surface 
field satisfies the criterion for the magnetic buoyancy instability, the field emerges again into the 
upper atmosphere after t ~ 4.1 hr. Eventually, several magnetic domes attain a height of ~ 7000 km 
at t ~ 5.6 hr. The size of the active region grows to 48 Mm x 12 Mm, and the photospheric flux 
amounts to 3.3 x 10 20 Mx. This two-step feature is consistent with our 2D models (papers I and II). 

We also observed multiple separation events of the magnetic elements at the photosphere, pos- 
terior to the horizontal escaping flow of the compressed plasma. Such separating elements move 
apart from each other at the rate of 4-8 km s" 1 and stop at the edges of the active region, which is 
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determined by the extension of the subpho to spheric field. The multiple separations are the results 
of the interchange-mode instability of the rising tube while in the convection zone. The magnetic 
elements then gather at the edges to make two alignments of the "pores" (tadpole-heads). The 
alignments show shearing motions at the rate of ~ 0.5 km s" 1 , which is explained by the inner 
field emergence (kinematic effect) and the Lorentz force effect (dynamic effect). Upflows of a few 
km s" 1 and downflows up to 1 km s" 1 are observed in the emergent areas and in the cores of sur- 
face fields respectively. As far as we know, we have never observed such photospheric features, 
especially the multiple separations, in previous 3D calculations applying flux tube as an initial con- 
dition (e.g. FanlboOl ; Archontis et alJboCM ). Note that some calculations using flux shee ts showed 
multiple separations via similar instabilities ( Isobe et al]|2005 ; Archontis & HoodlboOg K 



These f eatures are strongly rem iniscent of the observations of NOAA AR 5617 by 



Strous etal 



(1 1996b and 



Strous & Zwaanl (11999b . They found that faculae of both polarities separate from each 



other toward the edges of the region, which is similar to our findings of multiple separations. Also 
the shearing motions of the aligned tadpole-heads are consistent with the pores moving along the 
edges of the region toward the main sunspots. The difference of the tilt of magnetic elements 
between these two cases is caused by the twist direction. We assumed a right-handed tube in the 
initial state, which is favorable for the southern hemisphere. AR 5617 appeared in the northern 
hemisphere, which yields the left-handed twist. 

Th e su mmary of the comparis on between our results and the observations by 
d 19961) and 



Strous et al 



Strous & Zwaan 



(11999b is presented in Table \T\ The size of AR and the velocities are 



consistent between the two. As for the difference of separation speeds, note that lHarvey & Martin 
(11973b observed that faculae of opposite polarities separate initially at the rate of > 2 km s" 1 . The 
photospheric total flux is one digit smaller than the observed value, and we did not find any ma- 
jor sunspots in our active region. These differences may be because, in our calculation, the rising 
tube stops in 1.4 hours after it appears at the surface. The age of AR 5617 was estimated to be 
6.5-7.5 hours old at the beginning of the observation and 8-9 hours old at the end. In the obser- 
vation, emergence events showed undulatory structures with a typical wavelength of 8 Mm, while, 
in our case, each field line does not show undulation at around the surface. Such undulating fea- 
tures might be found in the more resolved calcu l ations or in the c alculations including thermal 



convection (llsobe et al 



2007 



Cheung et al 



2010). 



Strous & Zwaanl d 1999b summarized their ob- 



servations as a model in which each emergence event (separation) occurs in a single vertical sheet, 
which forms a series of sheets aligned in a parallel fashion (see Fig. 8 of their paper). From our 
results, the separations are explained as the consequence of the interchange-mode instability of the 
flattened flux tube beneath the surface (see e.g. Fig.0. 

On the basis of the numerical results in this paper, we suggest a theoretical picture of 
the flux emergence and the formation of an active region through the surface, which includes 



Strous & Zwaanl dl999)'s model. Our model is schematically illustrated in Fig. [3(1) The flux 



tube rises through the convection zone due to magnetic buoyancy. On approaching the surface, the 
tube decelerates and becomes flattened because of the photosphere in front of the rising tube. (2) 
When the photospheric field satisfies the criterion for the magnetic buoyancy instability, further 
evolution to the corona breaks out. Due to the interchange-mode instability, photospheric magne- 
togram shows multiple separation events. The separated elements reach the edges of the region, 
whose size is determined by the subsurface field. These elements then make two alignments of 
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Table 1. Summary of Comparison with AR 5617 



Simulation Results 



AR 5617° 



size of the region 

vertical unsigned total flux 

age from the appearance 

separation speed 

shearing speed 

upflow velocity 

downflow velocity 

wavelength of emergence pattern 



48 Mm X 12 Mm 
3.3 x 10 20 Mx 
0-1.4 hr 
4-8 km s" 1 
~ 0.5 km s" 1 
a few km s _1 
< 1 krns" 1 



50 Mm x 30 Mm 
4 x 10 21 Mx 

6.5-9 hr 
0.84 km s" 1 b 
0.73 km s" 1 c 
0.86 krns" 1 d 
1.26 krns" 1 e 
~ 8 Mm 



Notes. (a) IStrous et al 



( 19961) and lStrous & Zwaanl dl999h . Separation of facular elements. (c) Separation 



of pores. (J) Upflow in an emergent region. (e) Downflow in a facula. 




Fig. 7. Picture of flux tube emergence and the formation of an active region derived from the cal- 
culation results in this paper. Solar surface is indicated by a horizontal slice. (1) Flux tube rising 
through the convection zone decelerates to form a flat structure beneath the surface. (2) Magnetic 
elements of opposite polarities are observed to separate due to the interchange-mode instability. 
White and black ellipses indicate the positive and negative elements. These elements gather at 
the edges of the active region to make two alignments of pores. (3) Aligned pores show shearing 
motion when inner fields emerge. 



pores. Compare this panel with Fig. 8 of IStrous & Zwaanl (1 19991) . (3) As the emergence continues, 
inner fields of the flux tube rise and footpoints shift to show shearing motions. Lorentz force also 
drives the pores to shear. 
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The initial condition of the present experiment (# t ube = 67#o = 2.0 x 10 4 G, 7? tu b e 



5Ho 



1000 km, and q = 0.1 /Ho = 5.0 x 10 4 km 1 ) is t he same as that of Case 5 of our 2D cross- 
sectional calculation in lToriumi & Yokoyamal (1201 II) . One of the basic differences between these 



two simulation results is the rising time from the initial depth of -20 Mm to the surface. In 2D 
calculation, the flux tube reached the surface in t ~ 550ro = 3.8 hr, while, in 3D case, it took 
t ~ 600to = 4.2 hr. That is, the 3D tube rises slower. In 3D case, plasma in the tube apex drains 
down along the field lines to the both feet of the tube, which drives the tube more buoyant. At the 
same time, in the 3D regime, the magnetic curvature force pulls down the rising tube. The time 
difference between the two cases indicates that the curvature force dominates the draining effect in 
the present 3D experiment. 

In the solar interior as well as in the photosphere, several classes of convection may affect the 
rise of magnetic flux. Recent observ ations by SO T on board Hinode satellite have revealed the 
convective nature of the surface field ( Litesl2009 ). Stronger pores (> lkG) or sunspots, which are 
not found in this calcu lation without convection, would be formed through the convective collapse 
process JparkeJll978L Larg e-scale u pflow supports the rising of the tube, and more flux could be 



transported to the surface (IFan et al .112003b . Surface convection cr eates undulating fields, and the 
cancellation of such fields removes the mass trapped in the U-loop (Isobeetal ]2007ilcheung et al 



20101) . This process accelerates the density draining from the surface layer, which may be important 



Stein et al 



(12011 ) reported that flux of 20 kG in their 



to the spot formation. On the other hand, 
convective experiment, the same as used in ours, is too strong, because it produces large, hot, bright 
granules at the surface, which are not seen in the Sun. This is an interesting difference between the 
two types of calculations. It should be noted that the initial settings are different in th e two cases; w e 



, in lstein et 



al 



(2011)'s 



used an initial horizontal flux tube in the convectively stable interior, while, i 
calculation, uniform horizontal fields are advected into the computational domain by the convective 
upflows from the bottom boundary. More theoretical and observational studies are needed on the 
effects of the convection in the whole flux emergence process. 

Acknowledgements. Numerical computations were carried out on Cray XT4 at the Center for Computational Astrophysics, 
CfCA, of the National Astronomical Observatory of Japan. The page charge of this paper is supported by CfCA. This work 
was supported by the JSPS Institutional Program for Young Researcher Overseas Visits, and by the Grant-in- Aid for JSPS 
Fellows. We thank Dr. Y. Fan of the High Altitude Observatory, the National Center for Atmospheric Research, and the 
anonymous referee for improving this paper. We are grateful to the GCOE program instructors of the University of Tokyo 
for proofreading/editing assistance. 



References 

Archontis, V. & Hood, A. W. 2009, A&A, 508, 1469 

Archontis, V., Moreno-Insertis, F., Galsgaard, K., Hood, A., & O'Shea, E. 2004, A&A, 426, 1047 
Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford Clarendon Press, London) 
Cheung, M. C. M., Rempel, M., Title, A. M., & Schiissler, M. 2010, ApJ, 720, 233 
Fan, Y. 2001, ApJ, 554, LI 11 

Fan, Y, Abbett, W. P., & Fisher, G. H. 2003, ApJ, 582, 1206 

Harvey, K. L. & Martin, S. F. 1973, Sol. Phys., 32, 389 

Isobe, H., Miyagoshi, T., Shibata, K., & Yokoyama, T. 2005, Nature, 434, 478 

Isobe, H., Tripathi, D., & Archontis, V. 2007, ApJ, 657, L53 

Lites, B. W. 2009, Space Sci. Rev., 144, 197 

Manchester, IV, W. 2001, ApJ, 547, 503 



11 



S. Toriumi & T. Yokoyama: 3D MHD simulation on the solar flux emergence 

Matsumoto, R., Tajima, T., Shibata, K., & Kaisig, M. 1993, ApJ, 414, 357 
Newcomb, W. A. 1961, Physics of Fluids, 4, 391 
Parker, E. N. 1955, ApJ, 121, 491 
Parker, E. N. 1966, ApJ, 145, 811 
Parker, E. N. 1978, ApJ, 221, 368 

Shibata, K., Tajima, T., Steinolfson, R. S., & Matsumoto, R. 1989, ApJ, 345, 584 

Stein, R. E, Lagerfjard, A., Nordlund, A., & Georgobiani, D. 2011, Sol. Phys., 268, 271 

Strous, L. H., Scharmer, G., Tarbell, T. D., Title, A. M., & Zwaan, C. 1996, A&A, 306, 947 

Strous, L. H. & Zwaan, C. 1999, ApJ, 527, 435 

Toriumi, S. & Yokoyama, T. 2010, ApJ, 714, 505 

Toriumi, S. & Yokoyama, T. 2011, ApJ, 735, 126 

Zwaan, C. 1985, Sol. Phys., 100, 397 



12 



